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We consider a series of distorted black hole initial data sets, and develop techniques to evolve them 
using the linearized equations of motion for the gravitational wave perturbations on a Schwarzschild 
background. We apply this to 2D and 3D distorted black hole spacetimes. In 2D, waveforms for 
different modes of the radiation are presented, comparing full nonlinear evolutions for different 
axisymmetric i— modes with perturbative evolutions. We show how axisymmetric black hole codes 
solving the full, nonlinear Einstein equations are capable of very accurate evolutions, and also how 
these techniques aid in studying nonlinear effects. In 3D we show how the initial data for the 
perturbation equations can be computed, and we compare with analytic solutions given from a 
perturbative expansion of the initial value problem. In addition to exploring the physics of these 
distorted black hole data sets, in particular allowing an exploration of linear, nonlinear, and mode 
mixing effects, this approach provides an important testbed for any fully nonlinear numerical code 
designed to evolve black hole spacetimes in 2D or 3D. 
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I. INTRODUCTION 



As numerical relativity is empowered by ever larger computers, numerical evolutions of black hole data sets are 
becoming more and more common The need for such simulations is great, especially as gravitational wave 

observatories like the LIGO/VIRGO/TAMA/GEO network are gearing up to collect gravitational wave data over the 
next decade (see, e.g., Ref. Q and references therein). A recent, thorough, and very detailed study indicates that 
black hole collisions are considered a most likely source of signals to be detected by these observatories . Among 
the conclusions of this work is that reliable information about the merger process can be crucial not only to the 
interpretation of such observations, but also could greatly enhance the detection rate. Therefore, it is crucial to have 
a detailed theoretical understanding of the coalescence process, and particularly the merger phase, that can only be 
achieved through numerical simulation. 

However, numerical simulations of black holes have proved very difficult, largely because of the problems associated 
with dealing with singularities present inside. Even in axisymmetry, at present it is difficult to evolve black hole 
systems beyond about t — 150Af, where M is the mass of the system In 3D calculations, the huge memory 
requirements make these problems much more severe. The most advanced 3D calculations based on traditional Cauchy 
evolution methods published to date, utilizing massively parallel computers, have difficulty evolving Schwarzschild 
Misner or distorted Schwarzschild beyond about t = 50Af. Characteristic evolution methods have been used to 
evolve distorted black holes in 3D indefinitely although to date waveforms have not been extracted and verified, 
and it is not clear whether the technique will be able to handle highly distorted or colliding black holes due to potential 
trouble with caustics. 

In spite of these problems, much physics has been learned, especially in axisymmetry. There, calculations of distorted 
black holes with |£-11| and without angular momentum [ p2| , p^, h ave been performed. Furthermore, calculations of 
the Misner two black hole initial data have been carried out pDlJ], and the waveforms generated during the collision 
process have been extensively compared to calculations performed using perturbation theory p"5|-p7[. One of the 
important results to emerge from these studies is that the numerical and perturbative results agree very well, giving 
great confidence in both approaches. In particular, the perturbative approach turned out to work extremely well in 
some regimes where it was not, a priori, expected to be accurate. This has led to considerable interest in comparing 
perturbative calculations to full scale numerical simulations, as a way of not only confirming the validity of each 
approach, but also as an aid to interpreting the physics of the systems []l5|-|l7|. We expect this rich interplay between 
perturbation theory and numerical simulation to play an important role in the future of numerical relativity. 
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To this end, in the present paper we study a family of distorted axisymmetric and full 3D black hole initial data sets 
using perturbative techniques. These data sets consist of single black holes that have been distorted by the addition 
of gravitational waves. They can be considered to represent the final stages of the coalescence of two black holes, just 
after the horizons have merged. They therefore provide an excellent system to study the late stages of this process 
without having to model the long and difficult inspiral period. 

By using perturbative techniques to compute the waveforms expected in the evolution of these data sets, we provide 
an important testbed for full nonlinear codes that should evolve the same systems. In regimes where the distortions 
are considered moderately small perturbations of the underlying Schwarzschild or Kerr geometries, full scale numerical 
calculations should agree well with the type of perturbative treatment presented here. We expect that these results 
should have many uses in testing numerical relativity codes, and should also be useful in interpreting the physics 
contained in such simulations. For cases where different 3D codes and techniques are being developed, these results 
should help certify and interpret the numerical results. For example, 3D codes are being developed with standard 
ADM and new hyperbolic formulations of the equations [Q, with different slicing conditions. New approaches to 
black hole evolution that promise to extend dramatically the accuracy and duration of such calculations are under 
development, such as apparent horizon boundary conditions (see, e.g., 19,2^] and references therein) or very recently 
characteristic evolution 21 2^. In all cases, these testbeds should be applicable and should provide important tests 
of the numerical simulation results. 



II. THE PERTURBATIVE APPROACH TO EVOLUTION 



In this section we provide the basic mathematical formalism for evolving distorted black holes as perturbative sys- 
tems. We want to investigate the possibility of treating the non-spherical multipole moments of a general (numerical) 
spacetime as a linear perturbation about its "background" spherical part. Roughly speaking, this approach should 
be a valid way to describe black hole spacetimes whose nonspherical departure from Schwarzschild is small. In prac- 
tice, such an approach has already shown itself to be spectacularly successful in the "close-limit" approximation that 
regards the Misner initial data for two colliding black holes as a nonspherical perturbation to Schwarzschild | p7| , p3[ . 
In this paper we apply similar ideas to the evolution of distorted single black hole spacetimes, providing detailed 
comparisons with fully nonlinear numerical evolutions in axisymmetry, and providing the framework for comparing 
fully 3D simulations with perturbative evolutions. This approach is also motivated by earlier work of Ref. where 
perturbative evolution was used as a check on numerical evolution, but not as a way to evolve Cauchy initial data. 

We assume the spacetime metric gap to be in some general coordinate system {t,ri,0,(f)). This metric comes either 
from an analytic solution to the Einstein equations, or a numerical simulation. At this point we do not distinguish 
between these cases. In this coordinate system we can always write the metric in the form 

sph , 1 

where g''^^ is the spherically symmetric i — Q term, in a decomposition of gap into tensor spherical harmonics, and 
hap contains the higher multipoles which describe the deviations from the spherically symmetric i — Q mode. The 
terms hafj satisfy, to 0{h?)^ the linear field equations, linearized about g^^p ■ Thus if the terms hap are small, in the 

sense that hap <^ g^p^ they can be well approximated by a solution to the linearized equations. 

Our general approach then, is to take Cauchy initial data which we expect to describe a system which is close 
to spherical symmetry. From this initial data, we use a decomposition into tensor spherical harmonics to construct 
a background metric g^p, and perturbations hap. The perturbations hap are then evolved using the linear field 
equations. 

The methods for extracting g^^'' and hap, and for evolving hap are now discussed in detail. 



A. Extraction 



We write the background metric in \ 
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The perturbations /iq./?, which describe the deviations from spherical symmetry, are expanded using Regge- Wheeler 
harmonics |2j] as 

(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
(8) 

- coi6Yirn.4)/ sin6' (9) 
-4'") sin - cot eY,ra,e - sin^ e)/2 (10) 

/i^^ = sin^ er,,^ 

+i?2G(^™) + sine cos 

- c^^"' sin 6'(y£„, 60 - cot6Yi,rn,4,)- (11) 

In the above, there are seven even-parity {Hq™'\ h[^™\ h\^™^\ i?2^™'', hf™'\ K^^''^\ G*^^™-') and three odd-parity 
^^(£m) ^ ^(£m) ^ ^(£m)-j variables, which are functions only of t and 77. In Eqs. (||) to (^l|) a summation over the modes 
i> 1, —£< m < i is understood. 

Given this formal expansion of the full metric gap, using the orthogonality of spherical harmonics we can extract 
the components of the background metric 5^^'' by appropriate integrations over the 2-sphere 
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(We note that a similar expression in Ref. [p^ contained an error). The ten variables describing the non-spherical 
contributions (for £ > 2) can be extracted by similar integrals, for example 

Hi'"'\t,r^) = ^ J g^,Y;,^dn. (15) 

A complete list of formulae is given in Appendix In the usual case where the full metric gap is given numerically, say 
after solving the Cauchy initial value problem of numerical relativity, these integrals over 2-spheres can be computed 
numerically. In certain cases, the initial data may be known analytically, in which case the integrals can also be 
computed analytically. In the sections below we shall encounter examples of both cases. In any case, by performing 
these integrals over a series of 2-spheres of different radii on a given time slice (say, the initial time slice), we obtain 
the spherical background metric coefficients N , A, R and for each nonspherical ^m— mode the ten metric perturbation 
functions, e.g., H^™''^ ■ 

These metric perturbation functions, coming directly from the metric, are gauge dependent; their values will 
depend on the particular gauge used, say, by the numerical code used to generate them. However, a gauge-invariant 
perturbation theory of spherical spacetimes has been developed by various authors ||2^-|27[|. As shown originally by 
Abrahams [^-^, such formalisms can be used in numerical relativity calculations to isolate the even- and odd- 
parity gauge-invariant gravitational wave functions. The basic idea is that while the metric perturbation quantities, 
such as , will transform under infinitesimal coordinate transformations x^^ ^ + Sx^ in the usual way, one 

can use this information to construct special quantities that are invariant under such gauge transformations. On a 



3 



Schwarzschild background these gauge-invariant quantities are found to obey the standard Regge- Wheeler and ZeriUi 
equations describing gravitational waves propagating on the spherical black hole background. We will follow this 
approach below. 

To compute these gauge-invariant perturbations on the general, time- dependent background , we could follow 
the work of Refs. and construct gauge-invariant functions from the extracted i > 2 multipole moments, and 

evolve these functions using the gauge-invariant linearized equations. In this way we could treat data which is given 
in a general, and even time dependent, spherical coordinate system (although we note that in the present formalism 
the shift terms gu must be treated as perturbations, i.e., they must be formally 0{h)). Here, we take a simpler 
approach, which has so far been suitable for our needs. We assume that the Cauchy initial data for g^p is, to 0{h), 
given on a hypersurface of constant Schwarzschild time. Wc then construct gauge invariant variables using Moncrief 's 
prescription [p5|, which can then be easily evolved using the ZeriUi or Regge- Wheeler wave equations. 

We first transform gaf3 to the areal radial coordinate, using 




and then 
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Note that dR/drj can be easily calculated using 
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Then, we calculate all the required multipole moments of this transformed metric, integrating the metric over 2- 
spheres as described above, using formulae detailed in the Appendix This provides the ten metric perturbation 
functions as described above. 

With this information, following Moncrief we construct an odd-parity function r), and an even-parity function 

Q^„i{t,r), which are invariant under first order coordinate transformations. These are defined by 
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Here we also note the definition 

r 

where A'l is the Schwarzschild mass of the background. We approximate this mass by examining the function m(r), 
given by 

m{r) = ^-(^l-^^^M + 0{h). (18) 

In practice, in the spacctimes we are considering here, the function m{r) is quite constant in r, and so there is no 
practical ambiguity in computing the mass. Alternatively, for very small perturbations, one could consider simply 
using the ADM mass of the spacetime for M, as is often done jl2j, but this will include contributions from the waves 
in the spacetime. Hence the mass defined by Eq. ( p^ ) provides a slight improvement in the treatment. 

While the individual metric quantities hap will transform in the usual way by a linear gauge transformation 
x'^ ^ x'^ + Sx^ , the functions Q^^_^ and are invariant, and thus are more directly connected to the physics of the 
system. 
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We note that we have introduced a sUght inconsistency in our construction of the gauge-invariant quantities 
and Qp^- We have defined the Regge- Wheeler perturbation functions (e.g. H^^"^^) in a general way, on a general 
time-dependent background metric, but in computing the gauge-invariant quantities we have simply assumed the 
background to be in Schwarzschild coordinates. More complex expressions for the gauge-invariants on the more 
general background could be used, but in the present applications this has not proved to be necessary. 

At this stage we hint at a complication in this procedure which will we return to in more detail below. The numerical 
prescription outlined above, using integrals over 2 spheres to pick off the different £ — m wave modes in a distorted 
black hole spacetime, simply lumps "everything not spherical" into the wavefunctions. It does not discriminate 
between contributions at various orders in the perturbation expansion, but rather it lumps them altogether. One 
must be carefiil when using the standard perturbation equations, as described below, which are derived from a formal 
theory keeping only terms at linear order. In particular, we will encounter cases where some perturbation modes 
appear only at higher order. 



B. Evolution 



Having taken a general distorted black hole metric, we have detailed a technique to compute the gauge-invariant 
functions describing the waves on a Schwarzschild background. We now turn to the linearized evolution of these 
quantities. These gauge invariant functions obey the wave equations 
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where the Regge- Wheeler potential VRwif) is given by 
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and the Zerilli potential is 
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and finally k = {I - l){l + 2) + QM/r and r* = r + 2Mln(r/2M - 1). 

In general, the first time derivatives of and must also be calculated to provide Cauchy data for the 
evolution. In this paper only time symmetric initial data has been used, for which 
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= 0, 



(23) 



but it is no problem to provide time derivatives through analysis of the extrinsic curvature variables given as initial 

data in the more general non-time symmetric case. 

To summarize the development of this section, we have detailed the technique we use to isolate the even- and 
odd-parity gauge-invariant perturbation functions and (5£,„ from a numerically generated metric, under the 
assumption that the waves are linear perturbations on a spherical background metric. This information can be used 
in two ways: (a) it can be used to extract waveforms from a numerical simulation at a finite radius, and (b) it can 
be used to provide initial data for the linearized evolution equations given originally by Regge- Wheeler (odd-parity) 
and Zerilli (even-parity). These evolution equations can be used to compare results obtained with the full nonlinear 
evolution of black hole spacetimes. 



III. APPLICATION TO DISTORTED BLACK HOLES 



In this section we take the somewhat abstract discussion of the previous section and show how one applies it in 
practice to an actual family of numerically generated distorted black hole data sets. The method is tested on the 
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so-called black hole plus Brill wave spacetimes |31|,^, which represent a black hole distorted by a gravitational wave. 
The deviation from a Schwarzschild spacetime can be parameterized by a dimensionless parameter a, corresponding 
to the amplitude of the wave. In the sections below we describe these initial data sets, describe tests we can perform 
on the procedures used to obtain initial data for the perturbation equations discussed above, and in Section IV we 
compare linear and nonlinear numerical evolutions. 

A. The Black Hole Initial Data 

In this section we review the single distorted black hole initial data sets that we evolve in this paper. These black 
holes are distorted by the presence of an adjustable torus of nonlinear gravitational waves. The amplitude and shape 
of the initial wave can be specified by hand, as described below, and a range of initial data can be created, from 
slightly perturbed to very highly distorted black holes. Such initial data sets, and their evolutions in axisymmetry, 
have been studied extensively, as described in Refs. Three dimensional versions have been developed and 



studied recently in Refs. [0,p2 33 



Following |31| ], we write the 3~metric in the form originally used by Brill | ]34| : 

d£^ = i}'^ (e^i [dif + dO'^) + sin^ Odcj)'^) , (24) 

where 77 is a radial coordinate defined by f = "t^''' '^here M and f are the mass and standard Schwarzschild isotropic 
radial coordinate of the black hole in the Schwarzschild limit, described below. In this paper, we choose our initial slice 
to be time symmetric, so that the extrinsic curvature vanishes, although more general datasets can also be considered 

ii- 

Given a choice for the "Brill wave" function g, the Hamiltonian constraint leads to an elliptic equation for the 
conformal factor tp. The function q represents the gravitational wave surrounding the black hole, and is chosen to be 

q{il,e,(l)) = asin'^e* (e"(^)' +e"(^)') {l + ccos^cl)) . (25) 

Thus, an initial data set is characterized by the parameters (a, 6, w, n, c), where, roughly speaking, a is the amplitude 
of the Brill wave, h is its radial location, w its width, and n and c control its angular structure. The parameters n, 
which must be a positive even integer, and c dictate the angular [6 — 0) structure of the Brill wave. In particular, if 
c = 0, the resulting spacetime is axisymmetric. 

A study of full 3D initial data and their evolutions are discussed elsewhere [ ^2p^j35[ . If the amplitude a vanishes, 
the undistorted Schwarzschild solution results, leading to 

/,(o) _ 



= V2Mcosh(^^j . (26) 

Note that the initial data defined by Eq. ( p^ ) and ( p5| ) have octant symmetry (i.e., equatorial plane symmetry and 
discrete symmetry in the four quadrants around the z— axis). Together with the time-symmetry condition, this implies 
that Qf^{t, r) = and Q^^{t, r) is real, and only non-zero for even t and m. Therefore in the tests performed in this 
paper, only even-parity modes with even ^,to numbers are present, and in all analysis that follows linear evolutions 
will be carried out with the Zerilli evolution equation. However, the techniques presented are quite general, and can 
be used on a wider class of initial data sets than considered here, containing the full range of modes. All the data 
sets considered now, and in the following sections have M = 2, w = 1 and & = 0, that is, loosely speaking, the Brill 
waves initially have unit width, and are centered on the throat of an M = 2 Schwarzschild black hole. 

The elliptic equation to be solved for tp must in general be solved numerically, as will normally be the case. Once 
this is done, it can be evolved with a fully nonlinear numerical relativity code, and it can be also be used to compute 
initial data for the perturbation equations as described above. 

However, in this case we can also write down a solution in terms of special functions which analytically solves the 
Hamiltonian constraint to linear order in the expansion parameter a. This will provide a useful check on our procedures 
to compute linear initial data for the perturbation equations from the full nonlinear, numerically generated data. We 
give these solutions in schematic form below, and in more detail in Appendix |b[ 

The conformal factor is expanded into spherical harmonics, such that 

00 I 

V;(77,0,0)=^(°)('7)+aE E ^2{v)Yi,n{0,cl>)+O{a^) (27) 

1=0 m=-i 
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and the Hamiltonian constraint is solved to linear order in the expansion parameter a. Depending on the parameters 
considered in Eq. (^5|), the resulting linear conformal factor has a different angular structure. Three cases are 
considered in the following: 

1. n = 4, c = 0: The only non-zero coefficients are V'oo^ ^20^ ^^'^ V'io''- 

2. n = 2, c = 0: In this axisymmetric case, the only non-zero coefficients in the expansion are 'ip^]^ and ip^Q . The 

^^40^ term that one might expect is missing. Note that this implies that at linear order the £ = 4 Zerilli function 
vanishes, although it will appear at order a^. 

3. n — 4:, c =/= 0: In this non-axisymmetric case, the coefficients V'oo^ '^20' ^ ^22^ V^io^ ^''■a ^^"^ V'4^12j ^'"^ 
non-zero. In this case the •ip^l'^ contribution does not exist, and hence the £ ~ m ^ A Zerilli function will vanish 
at linear order. 

With these fully nonlinear numerical solutions and perturbative expansions we are in a position to extract and 
analyze initial data for the Zerilli evolution equation (and for the full 2D nonlinear evolution code), which is the topic 
of the next section. 



B. Extracting Wave Initial Data 

In this section we discuss specific examples of the extraction procedure applied to the axisymmetric black hole 
initial data sets discussed above. Here we focus on the extraction of initial data itself, not on the evolution. We use 
the exact analytic solution to the perturbative initial data equations as an aid to evaluating and understanding the 
numerically extracted initial data from the full nonlinear solution to the Hamiltonian constraint. 

There are two main variations on the extraction of initial data that we compare and contrast in this section: 

(a) We compute the gauge-invariant functions numerically, by computing numerical integrals over 2-spheres 
according to the procedure outlined above and in Appendix ^ 

(b) We can use the exact perturbative solution to calculate analytically the gauge-invariant variable denoted by 
Qtm^ ^ which we obtain by linearizing about the spherically symmetric background, and then constructing the variables 
as if the background was Schwarzschild. This follows the spirit of the general numerical procedure described above, 
and is the solution to which we expect to converge with the complete numerical procedure. This exact solution is too 
lengthy to be reproduced here, but can be obtained straightforwardly with a computer algebra package. 

We note that we can also use the exact solution to the perturbation expansion to calculate a gauge-invariant variable, 
which results from linearizing the spacetime about the background mass M Schwarzschild spacetime. This is not quite 
the same as the wave Q^^^ obtained by linearizing about the more general spherical background. This is due to 
two effects. First, the individual metric perturbation functions, such as defined via Eq. (|l5|) above, are computed 
on a more different spherical background than Schwarzschild Second, the Brill wave introduces the presence of a 
non-zero coefficients ij^QQ in the expansion (|2^), which is reflected in the definition of the mass defined in Eq. (p^). 
The difference between the two waveforms is O(a^). 



1. Testing with Linear Solution 

The general method, and code, were tested using the exact linear solutions for initial data using the conformal 
factor described above. Using these exact linear solutions, the 3-metric components ( |2^ ) were constructed on a 2D 
grid in {rj,9) coordinates (keeping only those terms linear in the Brill wave parameter a). The waveforms QJq' Qto 
and Qqq were numerically extracted, as described in Sec. [I A, and compared to the analytically computed functions 



Q^rn^ ■ expected from the numerical methods used, the numerically extracted waveforms converged, as ©(AS^), 
to Q^^'* , for each £-mode considered. This provides an excellent test of the accuracy of the integration routines on 
the 2-sphere required to compute the gauge- invariant wavefunctions, and of the rather complex expressions involving 
the various tensor-spherical harmonics for different ^-modes that must be coded. 

As an example. Fig. ^ shows the numerically extracted and exact waves for the initial data sets (a) (a = 0.05, n = 
4, c = 0) and (b) (a = 0.05, n — 2^c — 0), plotted as a function of the logarithmic radius rj. The numerically and 
analytically computed wave, calculated from the perturbative analytic solution to the Hamiltonian constraint, are 
shown as the dotted and dashed lines in the figure. On this scale they are almost indistinguishable, demonstrating the 
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FIG. 1. Numerically extracted and exact initial data for the cases (a) n = 4, a = 0.05, c = 0, and (b) n — 2, a = 0.05, 
c = 0. In both graphs: the dashed lines show waveforms extracted exactly about the spherical background and the dotted lines 
waveforms extracted numerically from linear initial data. These are practically indistinguishable since they lie almost on top of 
each other. The solid lines are waveforms extracted numerically from nonlinear initial data. The solid vertical line indicates the 
position of the peak of the background Schwarzschild potential in this radial coordinate. All the waveforms have been scaled 
by the appropriate extracted mass. 

high accuracy of the extraction procedure. For comparison, we also show the numerically extracted wave, discussed 
fully in the next section, computed from the full nonlinear solution to the Hamiltonian constraint, as a solid line. 
This indicates the difference between the linear and nonhnearly computed initial data. 

In Fig. pi, we show the ^ = 2,4,6 modes. For this initial data set with n = 4, the ^ = 2,4 modes are present at 
linear order, while the ^ = 6 mode is not. The numerically computed mode from the linear solution should vanish, 
and does so at second order in AO. Correspondingly, for the n — 2 linear initial data set, only the I ~2 mode should 
be present, and this is seen in Fig. ^d, which shows the numerically extracted 1 = 2 mode to be indistinguishable 
from the exact result, and the £ = 4 mode to be very close to zero. 



2. Extractions from the Nonlinear Initial Data 



Now we use the procedure detailed in Sec. |l| to extract waveforms from (2D numerical) initial data generated 
by solving the nonlinear constraints, for the same Brill wave parameters as above in Sec. IIIB 1. In Fig. ^ the initial 
waveform obtained is compared with the previous waveforms from the exact linear initial data. The waveforms match 
very closely, with the greatest differences near the black hole throat. This area, with the greatest difference between 
the waveforms from the linear and nonlinear initial data, is also inside the maximum of the background black hole 
potential, which is located (for ^ = 2) at ry sa 1.37, and most of this part of the waveform will be radiated across the 
horizon to disappear in the black hole. 

As the Brill wave parameter a is increased, the higher order terms begin to have more influence, and the waveform 
obtained from the nonlinear initial data deviates more from that from the linear initial data. This is demonstrated 
in Figs. 1^ and ^, in which as extracted from nonlinear data is shown for a range of Brill wave amplitudes, a. 
The Figures show that in the region close to the black hole throat (at r/ = 0) the amplitude of does not increase 
linearly with a, as expected. However, at the maximum of the black hole potential, r] « 1.37 (for 1 = 2 and M = 2), 
Qtm scales nearly linearly with a. 

We will see in the next section that even in cases such as these where the linear study deviates significantly from 
the nonlinear extraction, the linear evolution can still do a reasonably good job in predicting the results of the full 
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FIG, 2, Extracted waveforms from nonlinear initial data, compared to the corresponding exact linear solution for 
n = 4, c = for the different amplitudes a = 0,05, 0,1, 0,2, and 0,4, The solid lines show the waveforms extracted from the 
nonlinear initial data, and the dotted lines the exact linear waveforms. All the waveforms have been scaled by the appropriate 
extracted mass, and increasing size of a corresponds to increasing amplitude of the waveform. The solid vertical line indicates 
the position of the peak in the Schwarzschild potential. 



nonlinear evolution far from the black hole. The key reason for this seems to be that the nonlinearities of these 
particular datasets are largest well inside the peak of the potential barrier, and hence they are kept inside and 
propagate down the black hole. This point will be discussed further in the next section. 

In this section we have demonstrated that our extraction method is accurate by comparisons with known linear 
solutions, and shown that the nonlinear initial data sets, parameterized by the Brill wave parameter a, yield initial 
waveforms which agree closely with the known linear solution for small a. In the next section we turn to evolutions 
of these numerically extracted linear initial data sets. 



IV. EVOLUTION OF AXISYMMETRIC INITIAL DATA SETS 

The waveforms (9^„, extracted from data on the initial hypersurface, can be numerically evolved using the Zerilli 



equation (19), We stress that in this case, we compute from the nonlinear computed initial data, i.e. from 



the full solution to the Hamiltonian constraint. In the absence of an exact solution for the perturbative initial data 
expansion, as will normally be the case, this will be the only way to obtain initial data for the linearized evolution 
equations. For the cases considered in this paper, only time symmetric initial data was considered, with = on 
the initial hypersurface. For non-time symmetric data, can be calculated using the numerical extrinsic curvature, 
using a procedure similar to that for extracting Q^"^ from gij , 

In the following sections, we compare these ID linear evolutions, with the results from 2D fully nonlinear evolutions 
of the original initial data sets [Q^, The data is compared by constructing, in both simulations, Q^^it) at the 
same Schwarzschild radius. In the nonlinear code (5^(t) this requires the same extraction procedure that was used 
to find the linear initial data, to be applied at the chosen Schwarzschild radius throughout the evolution. In all the 
examples shown here, the radiation waveforms were extracted at a Schwarzschild radius r = 15M, 

In principle one must be careful in comparing waveforms measured with different time coordinates. Since the 
nonlinear evolution implements a general slicing condition (in these examples an algebraic slicing), and not the 
Schwarzschild slicing used in the linear evolution, it could happen that corrections would be needed to account for 
the differences in slicing. However, in our simulations the slicings are similar enough in the regions where the waves 
are computed that this is only a very small effect, as is borne out by the Figures, 
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FIG. 3. Extracted waveforms from nonlinear initial data, compared to the corresponding exact linear solution Q'l^^ for 
n = 2, c = for the different amplitudes a = 0.05, 0.1, 0.2, and 0.4. The solid lines show the waveforms extracted from the 
nonlinear initial data, and the dotted lines the exact linear waveforms. All the waveforms have been scaled by the appropriate 
extracted mass, and increasing size of a corresponds to increasing amplitude of the waveform. Note that in the bottom figure, 
the amplitude of Q'^q^^ is zero for all values of a. The solid vertical line indicates the position of the peak in the Schwarzschild 
potential. 




For all the results shown here, both for the linear and nonlinear evolutions, the accuracy of the results were carefully 
studied, to be sure that any differences between the nonlinear and linear waveforms is due only to the different initial 
data or the mode of evolution. This is particularly important to ensure that the observed modes are not due to 
insufficient resolution or boundary effects. 



A. Evolutions for n=4, c=0 

First, we discuss the evolutions corresponding to the axisymmetric initial data sets specified by ti = 4, c = 0, for a 



range of Brill wave amplitudes a. For the exact linear initial data described in Sec. Ill A, it was seen that QJq and 
occurred to linear order in the Brill wave amplitude a, with all other modes occurring at higher order. This linear 
scaling of Q^q and Q^q was also seen in Fig. |2[ in the initial data used for these linear evolutions. 

Fig. ^ compares the radiation waveforms at r — 15M from the linear evolution of the Zerilli equation and from the 
nonlinear evolution of the full Einstein field equations. It is clearly seen that for the lowest amplitude Brill waves, 
the waveforms from linear and nonlinear evolutions match very closely. As the amplitude increases, the deviations 
become more apparent. 

In Fig. ^ we isolate the cases a = 0.05 and a — 0.4 for closer study, and compare the relative amplitudes of the 
£ — 2 and £ — A waveforms. This clearly shows the high level of agreement for the lowest amplitude, shown in Fig. ||. 
Note that the scales have been chosen in the two graphs, such that if the waves were scaled linearly with a, they 
would be the same size in both (a) and (b) . Although we have evolved the same (non- linear) initial data in both cases, 
nonlinearities are clearly coming in, even in the initial data. 



B. Evolutions for n=2, c=0 



The second group of initial data sets were created using n = 2 and c — 0, for the same range of Brill wave amplitudes 
as in the previous section. For these parameters, the exact linear initial data sets of Sec. Ill A contained only the QJq 
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FIG. 4. Comparison of (. — 2 and £ = 4 waveforms a,t r — 15M from linear (dotted line) evolutions of linearized initial data 
and nonlinear (solid line) evolutions of nonlinear initial data. The initial data were constructed with the Brill wave parameters 
n = 4, c = 0, for four amplitudes, a = 0.05, 0.1, 0.2 and 0.4. 
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FIG. 5. Comparison of (. — 2 and £ = 4 waveforms at r = 15M from linear (dotted line) evolutions of linearized initial data 
and nonlinear (solid line) evolutions of nonlinear initial data. Figure (a) corresponds to the initial data set n = 4, c = and 
a = 0.05, and Figure (b) corresponds to the higher amplitude case n = 4, c = and a = 0.4. 
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mode to linear order in a, with all other modes occurring at second order or higher in a. This behavior is also seen 
in Fig. 0, in the initial data used for these linear evolutions. 

Fig. compares the radiation waveforms at r = 15M from the linear evolution of the Zerilli equation and from the 
nonlinear evolution of the full Einstein field equations. For these initial data, the £ = 2 waveforms shown in Fig. ^ 
again match well for low amplitudes, with deviations growing for the higher amplitudes. The i — A waveforms, shown 
in Fig. ^3 are of much lower amplitude than the i = 2 waveforms. As is apparent in the Figure, these waveforms from 
the linear and nonlinear evolutions do not match well, in amplitude or phase for any of the Brill wave parameters used. 
However, we should not expect them to. We have already seen in Sec. [II A, in these n = 2, c = initial data sets, 
there is no £ ^ A content at linear order in the expansion parameter a. The expansion of the initial data in powers 
of the amplitude a is directly related to the formal expansion used in the derivation of the Zerilli evolution equation. 
The equation is valid only for evolving initial data of order a, but not order a^, since all such terms were dropped in 
the derivation. In the case in point, the £ = 4 data are of order and hence this data will not be accurately evolved 
with the linear evolution equation. Basically this means that the £ = A piece, being of order a^, is too small to be 
treated by linear theory, and only the nonlinear code can pick this up. 

As shown by Gleiser et al, js^, the Zerilli equation can be generalized to higher order, in which case it has the same 
form as for the linear case, but now with nonlinear source terms that mix in different multipoles. In this case we could 
in principle evolve the £ = 4 data with such an equation, taking into account quadratic terms in the linear £ = 2 (and 
other) modes. In practice however, this would not be possible with the extraction method as used here, since errors 
of order are already introduced. In fact, it is worth stressing that the standard linear extraction process used to 
obtain the waveforms for the nonlinear code has been applied even for these nonlinear modes. The implications of 
such higher order effects on waveform extraction should be carefully considered in future investigations. 

The main point to be made here is that linear theory does not, and should not agree with the nonlinear code in 
this case, and one must be careful in applying linear theory as a testbed. On the other hand, these techniques open 
the door to careful and systematic studies of nonlinear physics of black holes, such as mode mixing, that could be 
pursued in the future. 

Fig. shows the £ — 2 and £ = 4 waveforms on the same graph, for the low amplitude a = 0.05 and higher amplitude 
a — 0.40 cases. Fig. ^ shows results from the a — 0.05 initial data, showing very good agreement between the linear 
and nonlinear waveforms. The £ — 4 waveforms are too small to be seen on this scale. Fig. |^ shows results from the 
higher amplitude a — 0.40 case. Here the discrepancy between the two methods can be seen in both the £ — 2 and 
£ — 4 waveforms. 

One might expect, given the rather large discrepancy between linear and nonlinear initial data for the a = 0.40 
case, that the linear and nonlinear evolved waveforms would be rather more different than they are. However, as 
previously noted the largest deviations between the linear and nonlinear initial data occurs near the horizon, well 
inside the peak of the potential barrier. Although linear theory breaks down inside the potential barrier, linear theory 
is still fairly accurate outside t he p eak where waveforms are measured. This effect has also been seen in studies of 
black hole collisions (see, e.g. [^^-16; 2^,^) where perturbative treatment was very successful even when one might 
naively think it would break down. 



C. Radiated Energy 



The asymptotic radiated energy flux, for each angular mode, can be defined by pi 



dE- 



dt 




dQ 



dt 



dQ 



im 



dt 



(28) 



which can easily be numerically integrated to give the radiated energy, E''™ in each mode. Note that this energy 
formula comes from linear theory, and hence will not provide an accurate prediction for the energy in the higher order 
modes. 

Fig. 1^ uses log plots to compare the energies in the waveforms from the nonlinear and linear evolutions of the two 
initial data sets described above. For the n — 4 initial data sets, shown in Fig. the energy radiated in both £ ^ 2 
and £ = 4 waveforms is similar, for the two types of evolution and for the range of Brill wave amplitudes considered. 
The energy carried by the £ = 4 waveform is some 10 to 20 times smaller than the energy carried by the £ — 2 
waveform. The radiated energy scales as ©(a^), for both modes. 

In Fig. ^ we show the n = 2 family of initial data sets, and the picture is different, as was already seen in the 
evolutions. Whereas the radiated energy is similar from both methods of evolution for the £ = 2 waveform, the £ = 4 
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FIG. 6. Comparison of = 2 and £ — 4 waveforms at r — 15M from (a) linear (dotted line) evolutions of linearized initial 
data and (b) nonlinear (solid line) evolutions of nonlinear initial data. The initial data was constructed with the Brill wave 
parameters n = 2, c = 0, for four amplitudes, a = 0.05, 0.1, 0.2 and 0.4. 
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FIG. 7. Comparison oi £ — 2 and £ — 4 waveforms at r=15M from linear (dotted line) evolutions of linearized initial data 
and nonlinear (solid line) evolutions of nonlinear initial data. Figure a) corresponds to the initial data set n = 2, c = and 
a = 0.05, and Figure b) corresponds to the higher amplitude case n = 2, c = and a = 0.4. 
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FIG. 8. The radiated energy, calculated at r = 15M, using the radiation waveforms from both the linear (dashed lines) 
and nonlinear (solid lines) evolutions. The top graph is for the data set n = 4, c — 0, and the bottom graph is for the data set 
n = 2, c — 0. In both cases, the range of Brill wave amplitudes a = 0.05 to a = 0.40 was used. 




waveform is very different. The £ — A waveform resulting from the nonlinear evolution, contains a factor of around 60 
times more energy than the corresponding waveform from the linear evolution. The "radiated energy" in the £ = 4 
waveform scales as 0{a'^), for both methods of evolution. But we emphasize that we have used a linearized energy 
measure in this case, where nonlinear effects should be accounted for. 



V. 3D CALCULATIONS 



So far in this paper we have concentrated on axisymmetric initial data sets, in order to aid the understanding and 
interpretation of the method. In this section, the same procedure is applied to construct linearized initial data from 
a non-axisymmetric 3D initial data set. Our aim here is to show that the extraction and linear analysis techniques 
developed here carry over easily to the full 3D case, even though the extraction process is more complicated. The 
evolutions of these data sets, and comparisons with nonlinear calculations, are left for the next paper in this series [ p9[ . 

The 3D initial data set studied here again corresponds to a black hole distorted by a gravitational wave, belonging 
to the same family as described in Sec. [II A. The non-axisymmetry follows from choosing a non-zero parameter c in 
Eq. |2^. The nonlinear data was created using the parameters a = —0.1, c = 0.5, n = A, b = and w = 1. 

Again, an exact linear solution for the 3-metric was constructed for these parameters, as detailed in Appendix 
Starting with this exact 3-metric on the initial hypersurface, exact expressions for the waveforms Qf^''' were found, 
using the extraction procedure of Sec. ||. Using the exact linear solutions, the 3-metric was calculated numerically, on 
a 3D polar grid, in {ri,9,(j)) coordinates (keeping only those terms linear in the Brill wave parameter a). Waveforms 
for different £ and m were then numerically extracted and compared to the exact Qg^'^- As before, the numerically 
calculated waveforms converge to the exact waveforms as 0{A6^, AcfP) as expected from the numerical methods. 

Fig. ^ shows the numerically extracted and exact waveforms, for the modes which are present to 0{a) in the exact 
solution. In this figure, the top graph shows the higher amplitude axisymmetric (m — 0) modes, and the bottom 
graph the lower amplitude non-axisymmetric (m ^ 0) modes. The two solutions, shown by the dashed and dotted 
lines, cannot be distinguished for this resolution, {Nrj = 200,Ng = 52, — 52). 

Finally, instead of extracting the waveforms from the numerical linear initial data, we extract from 3D numerical 
data, found by solving the nonlinear constraints for the given Brill wave parameters. This 'nonlinear' initial data is 
constructed in the (rj, 9, 0) coordinate system. The waveforms obtained from this initial data are shown as solid lines 
in Fig. ||. Note that there are significant differences between the waveforms extracted from the linear and initial data. 
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FIG. 9. Numerically extracted and exact initial data for the case n = 4, a = —0.1, c = 0.5. The top and bottom show 
respectively the m = and m ^ waveforms which appear in the exact linear solution. In both graphs: the short dashed lines 
show waveforms extracted exactly about the spherical background; the dotted line waveforms extracted numerically from linear 
initial data; and solid lines waveforms extracted numerically from nonlinear initial data. All the waveforms have been scaled 
by the appropriate extracted mass, and the solid vertical lines indicate the location of the peak in the Schwarzschild potential 
barrier. 




However, the largest differences are again inside the maximum of the potential barrier of the background spacetime 
at ?7 « 1.37 for £ = 2. On the potential barrier the linear and nonlinear waveforms are very similar, with the largest 
difference occurring for £ = 2,m — 2, indicating that in an evolution of this time-symmetric initial data, there should 
be only small differences outside of this radius. Although not shown here, the difference between the initial waveforms 
from the linear and nonlinear data decreases as the Brill wave parameter a is reduced. 



Fig. 10 shows the initial waveforms which are present in the extraction from nonlinear data, but are not found in 
the linear initial data, where they occur at ©(a^) or higher. As in the axisymmetric cases discussed above, these 
modes should not agree with full nonlinear evolutions. 




0,5 1 1,5 2 2,5 

radial coordinalc, 77 

FIG. 10. Numerically extracted for the case n = 4, a = —0.1, c — 0.5. The graph shows waveforms extracted from the 
nonlinear initial data, for those £ and m modes which are not present in the exact linear solution. From top to bottom, the 
waveforms are Qji'j/M, Qf^/M, Q'^^/M, Qg^/M and Q^q/M. Here M is the extracted mass. 
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VI. CONCLUSIONS 



Building on previous work (e.g., [y_2[^3|3^) we have further developed the use of perturbation theory as a tool 
for numerical relativity. We presented techniques for studying single, distorted black holes without net angular 
momentum, using the perturbative Regge- Wheeler and Zerilli equations. These techniques include extracting initial 
data, decomposed into different i, m modes, from a numerically computed black hole initial data set, and evolving 
these modes forward in time with the perturbation equations. These results can then be used both to perform a 
powerful check on the accuracy of fully nonlinear evolution codes, from which waveforms can be extracted, and to 
help untangle linear and nonlinear effects, such as mode- mixing, of the full evolution. 

We applied this technique to a series of distorted black holes, using the "Brill wave plus black hole" family of initial 
data sets |3|,^. These data sets have adjustable parameters for the shape and amplitude of the initial distortions, 
and have been well studied previously and found to mimic the behavior of two black holes that have just collided 
head-on. They have a further nice property that an analytic solution to the initial value problem is possible if the 
Brill wave amplitude is small. This information was used both to test the accuracy of the perturbative initial data 
extraction techniques, and to understand the spectral decomposition of both the linear and nonlinear contributions 
to the initial data. 

We then performed a series of axisymmetric, fully nonlinear supercomputer evolutions of the full initial data sets, and 
compared gravitational radiation waveforms extracted from these simulations to results obtained via our perturbative 
techniques using the Zerilli equation. We find that for initial data sets containing £ ^ 2,4 modes that can be shown to 
be linear in the wave amplitude parameter, the extracted waveforms from the full nonlinear evolutions show excellent 
agreement with the perturbatively evolved waveforms, confirming both the very high degree of accuracy that can be 
achieved with axisymmetric codes and the linear perturbation code and approach to the problem. 

However, we also find that some modes computed with the nonlinear simulation code do not agree with the linearized 
treatment; in some cases the linear evolution can give very different results from the full nonlinear evolution for specific 
modes, even when other modes agree well. But in all such cases, one can see from analytic study of the initial data 
that it is precisely these "renegade" modes that show up at a nonlinear order in the wave amplitude parameter. Their 
perturbative evolution equations should then be extended to second order, with nonlinear source terms (composed 
of linear modes) coming in. In effect, such modes are too small to be treated by linear theory, since nonlinear 
contributions from other modes are relatively large enough not to be neglected. Hence one can see clearly the effects 
of mode-mode coupling in these cases. The combination of linearized treatment and perturbative treatment was 
essential not only to confirm the nonlinear code, but also to interpret the complex nature of the numerical results. 
Furthermore, we pointed out that the use of linear waveform extraction techniques should be studied further for 
systems that may contain modes only at higher, nonlinear order. 

Another important point to mention is the role that the potential barrier plays in the evolution of these distorted 
black holes. The application of perturbation theory to the evolution of these distorted black holes often works well, 
even in cases where one can see large deviation between linear and nonlinear treatment in the initial data. But in 
such cases, the error made by using a linear approximation is largest inside the peak of the potential barrier, and 
hence most of this nonlinearity is swallowed by the black hole. Hence, as has also been stressed in previous work , 
the potential barrier acts to trap much of the deviation from linear behavior, extending the range of applicability of 
perturbative treatment over what one might naively expect. 

Finally, we also applied this technique to study a new class of 3D distorted black hole initial data sets |Q , and we 
showed that the initial data can be extracted accurately for study in 3D, just as in axisymmetry. The comparison of 
evolutions of such 3D distorted black hole initial data between full nonlinear numerical relativity and perturbation 
theory will be the subject of the next paper in this series. 

In this paper we have only considered examples of time symmetric, even-parity perturbations of non-rotating black 
holes. The technique is much more general, and also applies to all even- and odd-parity modes, and to non-time 
symmetric initial data, which will be considered in future papers in this series. 

The data sets developed in Ref. js^] also contain distorted black hole initial data with angular momentum. The 
rotating case is considerably more complicated, and naturally involves using the Teukolsky equation to evolve pertur- 
bations of the Kerr metric. This important foUowup step will be considered in a future paper. 
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APPENDIX A: EXTRACTION FORMULAE 



Here we list the general formulae for extracting the Regge- Wheeler variables from the metric gai3 using the the 
spherical background g^^'' = di&g{-N{t, rf, A{t, rf, R{t, rf, R{t, rf sin^ 6), 



1 f I 

J ~ 9t0Y;^^^)dfl (Al) 



'+!) J sine' 



+— 




sm 




1 




e{i + 


1) 


1 




e(i + 


1) 



sin 



Hr'it,r) = — I g„Yl,dQ (A7) 



^i'"^ = 7(7^ / <i 9reYel,^e + -^^r^^-/™,^ } dn (A6) 

H['"'\t,r) = J g,,rY;^dn (A8) 

^^i'"^ (*. ^) = i2 / 5rrF/„rff2 (A9) 

- cot dQ (All) 

sni & J 

= !(i±l)G(^'»)(t,r) + ^1 (5^. + -^500) i^;™^^^^ (A12) 



APPENDIX B: LINEAR SOLUTION FOR CONFORMAL FACTOR 



In this appendix we derive an exact linear solution for the perturbed black hole initial data sets used in this paper. 
Expanding the conformal factor ip in terms of the Brill wave amplitude a, we write 

= + ai^(i) + 0(a2) (Bl) 

The zeroth order Hamiltonian constraint is then 

a^Vifo) + + cot edeip^°^ + csc^ 0dl^^°^ = (B2) 

with the spherical Schwarzschild solution 

= V2Mcosh(|) (B3) 
Keeping terms linear in a, the first order Hamiltonian constraint is 
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492i/;(i) +49^^(1) +4 cot 6iae7/;(i) +4 csc^ 619^1^(1) - = ~^'4j^°\dlq + dlq + 2 csc^ 6dlq) (B4) 
When ■(/j^^-' is expanded in spherical harmonics, 



(B5) 



=0 m=-£ 



Equation B4 reduces to a series ol ODEs, which can be solved for a particular case of the Brill wave function, q. The 
solution for two of the Brill wave functions considered in this paper are, 

1. g = 2asin^ 0exp(— 77^). 



The only nonzero coefficients in the expansion B5 are T/igg'' and which have the values 



^(1) ^ _ 2^/2M^ ^_^2 ^^gj^^^^2) + H!L[cosh(r;/2) + erf(7/) sinh(?7/2)] 



(1) 
20 



2VlOM7r 



2\/57r 
15 



e^/^ cosh(5V2) + \/5^15: 



,9/4 



erf ( -?7 + - I e 



-S-J/Z + erf ( 7? + - I e 



5»7/2 



(B6) 
(B7) 



2. g = 2asin''6le-'' (1 + ccos^ 0) 

For this case, which includes a non-axisymmetric contribution through the parameter c, there are several more 



modes in the expansion. The non-zero coefficients in B5 are now V'qq'', ^"20^ ^22 = ^^2^121 V'io^ ^^i2 = '^^^^-2- 



'A'20 



V'40 



— V27rM(2 + c) ^-cosh(^^j (2e-'^ + ^/¥ j + V^sinh (^^ j erf (r;) 



1 /27rM 



40. 



e-"' cosh (I) + 60Fe cosh ( y ) - U^/^e^/^ cosh ( ^ 



5?7 



+3\/^ei-^''/2erf (r; - 1) - 3V^ei+^'"/2erf (ry + 1) - 7^e^'^-'^''''^exi [ 77 - - 



7V¥e9/4+5''/2erf f 77 + I 
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cosh (I) - 44V^ecosh - U^e"^ cosh 



70 V 15 ^ 

-22V^ei-5''/2erf (r; - 1) + 22V^ei+5''/2erf (7/ + 1) - 7V^e^/4-5'7/2e,f j',^ _ 3\ 7^39/4+5^2^^^ j",^ + 3 



^ V2;^M (2 + c) 46-"' cosh (I) - 2V^e25/4 cosh (^^ j - ^e''>'^-'''^'^cri (^77 - ^ 

-p25/4+9r,/2g^f (^77 + ^ 

46""' cosh (I) - 2V^e25/4 cosh (^^^ - V^e25/4-9')/2erf (^77 - ^ 



1 /ttM 



If the parameter c = 0, the coefficients with ttt, 7^ vanish. 

Note that if the physical metric is required to linear order, (as used in the examples in this paper), then the 
components will be given by, for example 



= (1 + 2g)(7^("))4 + 4a(V;(^))' + ©(a^) 



(B13) 
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APPENDIX C: NUMERICAL METHODS 



The numerical implementation of the 2D and 3D nonlinear codes for the generation and evolution of initial data is 
described fully in §,|l|,|o|,|ll ■ 

The procedure for extracting waveforms in both the nonlinear initial data, and in the evolved nonlinear data is the 
same. The implementation consists of three main steps: 

(i) Find a 2-sphere on which the approximate spherical symmetry is manifest. For the examples used in this paper, 

appropriate spheres are known from the construction of the initial data, and the 2-spheres are simply spheres 
of constant isotropic coordinate radius rj from the center of the octant symmetry. 

(ii) Construct the spatial metric components (7^^, 5^0, .gr,^, gee, gecj,, gcji4, on the given 2-sphere. In general this 
procedure will involve interpolating metric components from the numerical grid used to create the initial data 
or for the evolutions to a 2D numerical grid on the given 2-sphere, and then transforming the metric to the 
polar coordinate system. However for the data sets used here, the initial data was created and evolved on the 
same 2-spheres used for the extraction. 

(iii) Integration over the 2-sphere to calculate the Regge- Wheeler functions, and directly from these the radiative 
variables and Qg^- The integrations were performed numerically using Simpson's rule, which calculates 
the integral with an error of C'(A0^,A0^). To avoid complications, the polar grid points on the sphere are 
arranged so that the polar axis is straddled. The number of polar grid points needed for an accurate estimate 
of the integral will depend in general on the angular structure of the metric components gaf} and the ^, m mode 
under consideration. 
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